## ----- REPLICATION CODE FOR ONLINE APPENDIX TO: CONTACT AND CONTEXT: HOW MUNICIPAL TRAFFIC STOPS SHAPE CITIZEN CHARACTER
## ----- AUTHORS: ALLISON ANOLL, DEREK EPP, AND MACKENZIE ISRAEL-TRUMMEL


## First run R Script from Article Replication 



## ----- Contents: ---------------------

#1. Cronbach's alpha on police evaluation variable for white and black subsets
#2. Alternative model specification (OLS)
#3. Zip code correlates
#4. Repeat models on racial subsets
#5. Interaction of investigative stops ratio and personal stop history
#6. Interaction of investigative stops ratio and proximal contact
#7. Alternative codings of investigative stops ratio: absolute value, dropping municipalities that are more likely to stop whites
#8. Subset to only those places with at least 10, 50, 100 stops of both whites and blacks. 


## ----- 1. Cronbach's alpha on racial subsets ---------------------

policeevaluation_whites <- cbind(ws$a005_1, ws$a005_2, ws$a005_3, ws$a005_4, ws$a005_5)
cronbach(policeevaluation_whites) ##Alpha  0.9068032

policeevaluation_blacks <- cbind(bs$a005_1, bs$a005_2, bs$a005_3, bs$a005_4, bs$a005_5)
cronbach(policeevaluation_blacks) ##Alpha 0.8985843



## ----- 2. Alternative model specification (OLS) ---------------------

evalindex_ols <- lm(index ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age + Illinois, data = stops)

vote_ols <- lm(voted ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + Illinois, data = stops)

nonvoting_ols <- lm(participation4 ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + Illinois, data = stops)


summary(evalindex_ols)
summary(vote_ols)
summary(nonvoting_ols)



## OLS regression table
texreg(list(evalindex_ols, vote_ols, nonvoting_ols),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age",
         "Illinois"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robustols",  stars = c(0.10, 0.05),   caption = "Policing Context and Personal Stop History on Views of Police and Participation, OLS", float.pos = "t", single.row = TRUE) 


## ----- #3. Zip code correlates  ---------------------

##adding zip correlates
zip1 <- read.csv(file="nhgis0012_ds233_20175_2017_zcta.csv")  # zip data from ACS, file 1
zip2 <- read.csv(file="nhgis0012_ds234_20175_2017_zcta.csv")  # zip data from ACS, file 2
names(zip1)
summary(zip1$ZCTA5A) #zip code
table(zip1$ZCTA5A)
zip1$zip <- zip1$ZCTA5A
summary(zip1$zip)

summary(zip1$AHY1E001) ##pop in zip
summary(zip1$AHZAE003) ##non-hispanic white in zip
summary(zip1$AHZAE004) ##non-hispanic black in zip
summary(zip1$AHZAE014) ##hispanic black in zip

zip1$white <- zip1$AHZAE003 / zip1$AHY1E001
summary(zip1$white)

zip1$black <- (zip1$AHZAE004 + zip1$AHZAE014) / zip1$AHY1E001
summary(zip1$black)

summary(zip1$AH04E001) ##total over 25 year olds in zip
summary(zip1$AH04E022) ##bachelor's degree in zip
summary(zip1$AH04E023) ##master's degree in zip
summary(zip1$AH04E024) ##professional degree in zip
summary(zip1$AH04E025) ##doctorate degree in zip

zip1$BA <- (zip1$AH04E022 + zip1$AH04E023 + zip1$AH04E024 + zip1$AH04E025)/zip1$AH04E001
summary(zip1$BA)

summary(zip1$AH1PE001) ##median hhi in zip in 2017 dollars
zip1$hhi <- zip1$AH1PE001
summary(zip1$hhi)

summary(zip2$ZCTA5A) #zip code
zip2$zip <- zip2$ZCTA5A
summary(zip2$zip)

summary(zip2$AH8VE001) #pop in zip
summary(zip2$AH8VE002) #us born in us
summary(zip2$AH8VE003) #us born in PR
summary(zip2$AH8VE004) #us born abroad

zip2$usborn <- (zip2$AH8VE002 + zip2$AH8VE003 + zip2$AH8VE004)/zip2$AH8VE001
summary(zip2$usborn)

summary(zip2$AIIFE001) #households in zip
summary(zip2$AIIFE002) #hh with cash assistance in zip

zip2$pubassist <- zip2$AIIFE002/zip2$AIIFE001
summary(zip2$pubassist)

zipdata <- data.frame(zip1$zip, zip1$black, zip1$white, zip1$hhi, zip1$BA, zip2$usborn, zip2$pubassist)
names(zipdata)
colnames(zipdata) <- c("zip", "zip.black", "zip.white", "zip.hhi", "zip.BA", "zip.usborn", "zip.pubassist")

stops.zip = merge(stops, zipdata, by="zip", all.x = T)
names(stops.zip)


##zip models
evalindex_zip <- lm(index ~ pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age + Illinois + zip.black + zip.white + zip.hhi + zip.BA + zip.usborn + zip.pubassist, data = stops.zip)

vote_zip <- lm(voted ~ pctstopinvest_ratio_zero + pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age + Illinois + zip.black + zip.white + zip.hhi + zip.BA + zip.usborn + zip.pubassist, data = stops.zip)

nonvoting_zip <- lm(participation4 ~  pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age + Illinois + zip.black + zip.white + zip.hhi + zip.BA + zip.usborn + zip.pubassist, data = stops.zip)

summary(evalindex_zip)
summary(vote_zip)
summary(nonvoting_zip)



texreg(list(evalindex_zip, vote_zip, nonvoting_zip),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age",
         "Illinois", 
         "Prop. Black", ##zip level variables
         "Prop. White",
         "Median Income",
         "Prop. BA or higher",
         "Prop. U.S.-born",
         "Prop. on Assistance"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:zipregs",  stars = c(0.10, 0.05),   caption = "Inclusion of Zip Correlates", float.pos = "t", single.row = TRUE) 


## ----- 4. Repeat models on racial subsets ----------------------


##black respondents
evalindex_bs <- lmer(index ~ pctstopinvest_ratio_zero + stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = bs)

vote_bs <- lmer(voted ~ pctstopinvest_ratio_zero + stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = bs)

nonvoting_bs <- lmer(participation4 ~  pctstopinvest_ratio_zero + stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = bs)

summary(evalindex_bs)
summary(vote_bs)
summary(nonvoting_bs)

##white respondents

evalindex_ws <- lmer(index ~ pctstopinvest_ratio_zero + stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = ws)

vote_ws <- lmer(voted ~ pctstopinvest_ratio_zero +  stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = ws)

nonvoting_ws <- lmer(participation4 ~  pctstopinvest_ratio_zero + stopped.bin + networkfel + victim + felcon  + income + ideo + education + woman  + age  + (1 | Illinois/zip), data = ws)

summary(evalindex_ws)
summary(vote_ws)
summary(nonvoting_ws)





texreg(list(evalindex_bs, vote_bs, nonvoting_bs),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped by police",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust3A",  stars = c(0.1, 0.05),   caption = "Restricted to Black respondents", float.pos = "t", single.row = TRUE) 



texreg(list(evalindex_ws, vote_ws, nonvoting_ws),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped by police",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust3B",  stars = c(0.10, 0.05),   caption = "Restricted to White respondents", float.pos = "t", single.row = TRUE) 




#----- 5. Interaction of investigative stops ratio and personal stop history ---------------

evalindex.int <- lmer(index ~ pctstopinvest_ratio_zero * as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

vote.int <- lmer(voted ~ pctstopinvest_ratio_zero * as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

nonvoting.int <- lmer(participation4 ~  pctstopinvest_ratio_zero * as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


summary(evalindex.int)
summary(vote.int)
summary(nonvoting.int)

## binary stops variable

evalindex.int.bin <- lmer(index ~ pctstopinvest_ratio_zero * stopped.bin  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


vote.int.bin <- lmer(voted ~ pctstopinvest_ratio_zero * stopped.bin  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


nonvoting.int.bin <- lmer(participation4 ~  pctstopinvest_ratio_zero * stopped.bin  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


summary(evalindex.int.bin)
summary(vote.int.bin)
summary(nonvoting.int.bin)


texreg(list(evalindex.int, vote.int, nonvoting.int),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age",
         "Stops Ratio * Stopped 1--2",
         "Stops Ratio * Stopped 3+"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust2A",  stars = c(0.1, 0.05),   caption = "Interacting Policing Context and Personal Stop History on Views of Police and Participation", float.pos = "t", single.row = TRUE) 




texreg(list(evalindex.int.bin, vote.int.bin, nonvoting.int.bin),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped by police",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age",
         "Stops Ratio * Stopped"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust2B",  stars = c(0.1, 0.05),   caption = "Interacting Policing Context and Binary Personal Stop History on Views of Police and Participation", float.pos = "t", single.row = TRUE) 





##----- 6. Interaction of investigative stops ratio and proximal contact ---------------

##create proximal contact binary variable
stops$proximal[stops$networkfel == 0] <- 0
stops$proximal[stops$networkfel > 0] <- 1
table(stops$networkfel, stops$proximal)


evalindex.intproxbin <- lmer(index ~ pctstopinvest_ratio_zero * proximal +  as.factor(stopped)  + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

vote.intproxbin <- lmer(voted ~ pctstopinvest_ratio_zero * proximal +  as.factor(stopped) + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

nonvoting.intproxbin <- lmer(participation4 ~  pctstopinvest_ratio_zero * proximal +  as.factor(stopped)  + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

summary(evalindex.intproxbin)
summary(vote.intproxbin)
summary(nonvoting.intproxbin)


texreg(list(evalindex.intproxbin, vote.intproxbin, nonvoting.intproxbin),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Proximal carceral contact",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age",
         "Stops Ratio * Proximal"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust2C",  stars = c(0.1, 0.05),   caption = "Interacting Policing Context and Proximal Contact", float.pos = "t", single.row = TRUE) 


#----- 7. Alternative codings of investigative stops ratio: ---------------

#create subset of only those places that are anti-black in their stop ratio

antiblackonly <- subset(stops, stops$pctstopinvest_ratio_zero > -0.000001)
summary(antiblackonly$pctstopinvest_ratio_zero)
length(antiblackonly$pctstopinvest_ratio_zero) #818 obs

#create an absolute value measure of bias (so anti-white and anti-black aren't differentiated)

stops$ratio_absval <- abs(stops$pctstopinvest_ratio_zero)
table(stops$pctstopinvest_ratio_zero, stops$ratio_absval)
summary(stops$ratio_absval)
length(stops$pctstopinvest_ratio_zero) #893 obs


#absolute value

evalindex.absval <- lmer(index ~ ratio_absval + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


vote.absval <- lmer(voted ~ ratio_absval + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)

nonvoting.absval <- lmer(participation4 ~  ratio_absval + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops)


summary(evalindex.absval)
summary(vote.absval)
summary(nonvoting.absval)



#anti-black only
evalindex.AB <- lmer(index ~ pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = antiblackonly)


vote.AB <- lmer(voted ~ pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = antiblackonly)


nonvoting.AB <- lmer(participation4 ~  pctstopinvest_ratio_zero + as.factor(stopped)  + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = antiblackonly)


summary(evalindex.AB)
summary(vote.AB)
summary(nonvoting.AB)


texreg(list(evalindex.absval, vote.absval, nonvoting.absval),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:",  stars = c(0.10, 0.05),   caption = "Investigative Stops Ratio Recoded as Racial Disparity", float.pos = "t", single.row = TRUE) 




texreg(list(evalindex.AB, vote.AB, nonvoting.AB),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:",  stars = c(0.10, 0.05),   caption = "Dropping Municipalities That Are More Likely to Stop White Motorists", float.pos = "t", single.row = TRUE) 




#----- 8. Subset to only those places with at least 10, 50, 100 stops of both whites and blacks  ---------------

##subset to only places with at least 10 stops of both whites and blacks in the data
stops_10blk <- subset(stops, stops$stopblk > 9)
stops_10 <- subset(stops_10blk, stops_10blk$stopwht > 9)
summary(stops_10$stopblk)
summary(stops_10$stopwht)


evalindex_10 <- lmer(index ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_10)

vote_10 <- lmer(voted ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_10)

nonvoting_10 <- lmer(participation4 ~  pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_10)

summary(evalindex_10)
summary(vote_10)
summary(nonvoting_10)


##subset to only places with at least 50 stops of both whites and blacks in the data

stops_50blk <- subset(stops, stops$stopblk > 49)
stops_50 <- subset(stops_50blk, stops_50blk$stopwht > 49)
summary(stops_50$stopblk)
summary(stops_50$stopwht)



evalindex_50 <- lmer(index ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_50)

vote_50 <- lmer(voted ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_50)

nonvoting_50 <- lmer(participation4 ~  pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_50)

summary(evalindex_50)
summary(vote_50)
summary(nonvoting_50)



##subset to only places with at least 100 stops of both whites and blacks in the data
stops_100blk <- subset(stops, stops$stopblk > 99)
stops_100 <-  subset(stops_100blk, stops_100blk$stopwht > 99)
summary(stops_100$stopblk)
summary(stops_100$stopwht)



evalindex_100 <- lmer(index ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_100)

vote_100 <- lmer(voted ~ pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_100)

nonvoting_100 <- lmer(participation4 ~  pctstopinvest_ratio_zero + as.factor(stopped) + networkfel + victim + felcon  + income + ideo + education + woman + white + age  + (1 | Illinois/zip), data = stops_100)

summary(evalindex_100)
summary(vote_100)
summary(nonvoting_100)


texreg(list(evalindex_10, vote_10, nonvoting_10),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust4A",  stars = c(0.10, 0.05),   caption = "Dropping Municipalities with $<$ 10 stops of Whites or Blacks", float.pos = "t", single.row = TRUE) 



texreg(list(evalindex_50, vote_50, nonvoting_50),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust4B",  stars = c(0.10, 0.05),   caption = "Dropping Municipalities with < 50 stops of Whites or Blacks", float.pos = "t", single.row = TRUE) 



texreg(list(evalindex_100, vote_100, nonvoting_100),
       custom.coef.names = c(
         "(Intercept)",
         "% Invest. Stops Ratio",
         "Stopped 1-2 times",
         "Stopped 3+ times",
         "Proximal carceral contact",
         "Crime victim",
         "Felony conviction",
         "Income",
         "Conservatism",
         "Education",
         "Woman",
         "White",
         "Age"
       ),
       custom.model.names = c("Evaluation of Police", "Turnout", "Participatory Acts"),
       dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:robust4C",  stars = c(0.10, 0.05),   caption = "Dropping Municipalities with < 100 stops of Whites or Blacks", float.pos = "t", single.row = TRUE) 


